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The formation of current sheets in ideal incompressible 
magnetohydrodynamic flows in two dimensions is studied 
numerically using the technique of adaptive mesh refine- 
ment. The growth of current density is in agreement with 
simple scaling assumptions. As expected, adaptive mesh 
refinement shows to be very efficient for studying singular 
structures compared to non-adaptive treatments. 



1. INTRODUCTION 

The formation of singularities in hydro- and magne- 
tohydrodynamic flows is still a controversial issue in the 
mathematics and physics community. Since mathemat- 
ically only very little is known Q, one has to rely on 
numerical simulations. Even in very elaborate numeri- 
cal experiments (see Bell and Marcus [Q], Kerr ||) non- 
adaptive treatment is limited very soon by the computer 
memory available, resulting in a resolution of less than 
512 grid points in each spatial direction. Since the sin- 
gular structures like tubes and sheets are not space fill- 
ing, adaptive mesh codes seem to be the right choice 
for studying these problems, as has been done by Pumir 
and Siggia |||| . Unfortunately, the methods used in [QJ|] 
could only refine the region around a singular point which 
lead in [Q to a substantial loss of energy. Of course, it 
is desired to refine all regions where the numerical reso- 
lution is insufficient. Modern adaptive mesh refinement 
algorithms, as introduced by Berger and Colella || and 
Bell et al. 0, do not possess the above limitations and 
are good candidates for studying singularity formation 
even in incompressible systems. 

In this paper, we investigate the formation of singular 
current sheets described by the ideal incompressible mag- 
netohydrodynamic equations (MHD equations) in two di- 
mensions for the time evolution of the velocity field u and 
magnetic field B. Using Elsasser variables z ± = u ± B, 
the MHD equations take the symmetric form 



d, 



Vz ± + Vp = , div z ± = 



(1) 



The equations (|l|) are integrated in a periodic quadratic 
box of length 2n using adaptive mesh refinement with 
rectangular grids self-adjusting to the flow. In each rect- 
angular grid, a projection method is used where the time- 
stepping is performed in a second order upwind man- 
ner [&]-p^]. For the projection step, we need the vortic- 
ities u) = (V x z ± ) • e z and potentials V' ± which are 
related by A?/^ = u> ± . 

The outline of the paper is as follows. In the next 
section, the adaptive mesh refinement algorithm is intro- 
duced. Then, we discuss the numerical results and com- 
pare the growth of current density with the prediction 
of Sulem et al. (llj . Finally, we conclude that adaptive 
mesh refinement is an ideal tool for studying singular 
structures and should be pursued further to study three 
dimensional problems as the finite time blow up in the 
incompressible Euler equations. 



2. ADAPTIVE MESH REFINEMENT 

2.1. General Strategy 

The main idea of adaptive mesh refinement is simple. 
One starts with a grid of given resolution and integrates 
the partial differential equation as usual. As soon as 
some criterion is fulfilled, this initial grid is refined. This 
is done by marking all critical grid points where the dis- 
cretization error exceeds a prescribed value. Then new 
grids with finer resolution and timestep are generated 
which cover all these critical points. These grids belong- 
ing to the next level arc then filled with interpolated data 
from the first level. Then, one integrates both levels un- 
til the resolution again becomes insufficient. Now the 
critical points are collected over all grids of the actual 
level being refined. Filling the new grids with data is 
achieved by first taking data from the previous level and, 
if existing, data from former grids of the same resolu- 
tion. This process is repeated recursively. In addition, to 
communicate the boundary conditions, each grid needs 
information about its parent grids and its neighbors. As 
one can see already, adaptive mesh refinement requires 
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the management of lists of levels, critical points, grids, 
parent grids and neighbors. Therefore, we programmed 
the handling of those structures in C++ whereas the nu- 
merically expensive integrations are done in fortran. In 
order to encourage the reader to use adaptive mesh re- 
finement we describe the above outline in more detail in 
the next paragraphs. 

To deal with all the different lists we defined tem- 
plated list classes and iterators which can be used for all 
classes representing levels, grids, critical points, parents 
and neighbors. 

The integrator used for all grids is based on a pro- 
jection method combined with second order upwinding. 
This scheme motivated by Bell et al. (|] was previously 
applied to incompressible magnetohydrodynamic flows in 
two dimensions fllO| ]. It is clear, that the equations un- 
der consideration can be easily exchanged by other ones 
using an explicit algorithm since the structures needed 
for adaptive mesh refinement and the integrator are in- 
dependent of each other. 

The timestep on a given level is advanced as illustrated 
by the following piece of pseudo-code. 

procedure integrate level 
do singlestep on level 
better boundary on level 
solve poisson equation on level 

if next level exists, then 

default boundary on next level 
do r times 

integrate next level 
update of level 

check criterion on level 

Starting at time to, the procedure singlestep performs 
one timestep Ati eve i on all grids of this level. In addition 
to the data within the grid itself the integration scheme 
needs boundary data which are by default obtained by in- 
terpolation in space and time from previous level data. In 
the subsequent procedure better boundary, the boundary 
data are, if possible, replaced by values from neighbor- 
ing fine grids. After boundary data have been commu- 
nicated on each grid, in order to perform the projection 
step Poisson equations with fixed boundary are solved 



for the potentials Atp 



± — / ,± 



Now all data of the actual 



level are advanced to a time to + Ati eve i and the recur- 
sion starts by integrating the next level if existing. The 
first step in this recursive process is achieved by supply- 
ing information about the default boundary data from 
parent grids. This is done by storing the increments cal- 
culated from the actual grids at time t and parent grids 
at time to + Ati eve \ . To achieve linear interpolation in 
time these increments are added to the boundary data 
at the end of singlestep. Storing only the increments in 
a special C++ boundary class avoids the memory over- 
head resulting from keeping data at present and previous 
times. On this next level, the timestep Ati eve i + \ and the 



spatial discretization lengths are divided by a refinement 
factor r. Therefore, the procedure singlestep has to be 
called r-times on this new level in order to reach the 
time to + Ati eve i . Having completed this recursive inte- 
gration loop, this level and all finer levels are advanced 
to time to + Ati eve i . Now the finer level data are used 
in procedure update to improve the values of the actual 
level. The procedure integrate is finished by checking if 
a certain criterion is fulfilled, that decides whether a re- 
finement step is performed. 



2.2. Regridding 

The criterion for refinement is adapted to the problem 
of current sheet formation. The global maximum of vor- 
ticity and current density is calculated and compared to 
the values when the last refinement was done. Regrid- 
ding is initiated, if the ratio of those maxima exceeds a 
prescribed value which is equal to the refinement factor r 
due to the scaling symmetry of the MHD equations (]]). 
The result of regridding is a new list of levels starting be- 
low the actual level. This new list replaces the old one, 
which is then deleted. 

The logical structure of the regridding procedure is 
shown in the subsequent pseudo-code. 

procedure regridding level 
for all grids on level 

mark critical points and append them to a list 
cover the critical points with rectangles (saw up) 

nesting rectangles into their parents and 
assign parents and neighbors 

fill the new rectangles with default data 

calculate global maxima for comparison 
in the procedure check 

if old level of same resolution existed before 

regridding, then 

better data on new level from old level 

solve poisson equation on new level 

if finer level existed before regridding, then 

regridding of new level 
else 

assign global maxima from old level 

else 

solve poisson equation on new level 

The procedure regridding starts with a loop over all 
grids of level to collect the critical points. Therefore, 
we calculate at each grid point the difference between 
the convection terms z T • Vz ± on the actual level and 
the next coarser one and if this difference exceeds a pre- 
scribed threshold e, we append this point and a surround- 
ing rectangle of given size to the list of critical points. In 
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the procedure saw up these critical points are covered 
with rectangles. The procedure nesting guarantees that 
they are properly nested into grids of the previous level 
allowing for more than one parent grid. At the same 
time, parent and neighbor grids are assigned to each new 
rectangle. Since these procedures are the most complex 
ones, they will be discussed in detail in the next subsec- 
tions. 

Now as each rectangle has information about his par- 
ents, the new rectangles are filled with spatially inter- 
polated data in default data. To avoid discontinuities, 
interpolation is done on the fields containing the highest 
derivatives, namely V xz ± . In order to supply boundary 
conditions for the solution of the Poisson equations, data 
for the potentials ip on the outermost boundary are as- 
signed as well. Afterwards, global maxima needed in the 
procedure check are calculated. 

If the recursive regridding was first invoked on the 
deepest level, the procedure is finished by solving the 
Poisson equations on the new level. Otherwise, data of 
the same resolution already existed and are used in better 
data to get more accurate values for the new grids. Data 
for the potentials are available after solving the Poisson 
equations. If the old level of same resolution was not 
the deepest level, the recursive regridding procedure is 
applied to the new level. In order to avoid unnecessary 
rebuilding of the level hierarchy, global maxima used as 
reference in check are assigned from the old level only in 
the other case. 



2.3. Grid Generation 

The grid generation is performed in the procedure 
saw up acting on a list of rectangles. On first entry, 
this list consists of one rectangle which covers all crit- 
ical points of that level. Each rectangle is now processed 
in the following way. First, it is decided in which direc- 
tion the first cut will take place. Therefore, we calcu- 
late vectors in the x- and y-direction which contain the 
number of critical points in each column or row, respec- 
tively. According to Bell et al. |7), we call them hori- 
zontal and vertical signatures E. The first cut is done in 
the direction with larger fluctuations in signature. This 
is achieved in the procedure cut dim, which first seeks 
for the best cut in this direction. In the procedure cut 
zeroes of the signature and its turning points (zeroes of 
A, = - 2E,- + E i+i) are taken into account as pos- 
sible cuts. If no such cuts are found, the mid point is 
chosen. A cut results in two lists of critical points. Each 
list is covered by a rectangle of minimal size. To ev- 
ery rectangle costs are assigned which are calculated as 
a sum of integration and memory costs (oc the area), 
boundary communication costs (oc the perimeter) and 
fixed costs (measuring the overhead for managing one 
additional grid). The two rectangles having the mini- 
mal costs are returned. Afterwards a loop over these two 



rectangles is performed. They are both given to the pro- 
cedure cut to find the best cut in the other direction. 
The costs of the two new rectangles in comparison to the 
original one's are used to decide whether the second cut 
is accepted or not. This gives a list of two, three or four 
rectangles. Their costs are summed up, and if they are 
less than the costs of the rectangle which entered the pro- 
cedure cut dim, they are returned to saw up. Otherwise, 
an empty list is given back. In the latter case, if the ef- 
ficiency measured by the ratio of critical points and grid 
points in the rectangle is insufficient, we enforce a cut in 
the middle of the longer side of the rectangle. Now the 
new rectangles are appended to a temporary list, which 
is, if not empty, passed to the recursive procedure saw 
up again. This recursion is stopped when further cuts 
do not allow a reduction of costs anymore. The above 
treatment is summarized in the two following pieces of 
pseudo-code. 

procedure saw up rectangles 
for all rectangles 

calculate E and variance in x- and y-direction 
if variance in x > variance in y, then 

apply cut dim on rectangle in ir-direction 
else 

apply cut dim on rectangle in y-direction 
if no cut found and efficiency insufficient, then 

half rectangle in longer direction 
append resulting rectangles to temporary list 
saw up of temporary list of rectangles 
if temporary list is not empty, then 

replace actual rectangle by temporary list 

procedure cut dim of rectangle in direction dim 
determine best cut in direction dim 
and return two rectangles 

loop over the two rectangles 
cut in other direction 

if costs are smaller than those of actual 
rectangle, then 

replace actual rectangle by list 

compare costs of new list (of 2-4 rectangles) with 
those of original rectangle and return cheapest 

An example, where saw up produces three new rectan- 
gles is shown in Figure |l]. 

2.4. Nesting 

After generation of non-overlapping rectangles in the 
procedure saw up, it is not guaranteed that all rectangles 
are properly nested in the rectangles of the parent level. 
A typical example, where this is not the case, is shown 
in Figure 2. 
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cuts of several levels. Therefore, we seek for cuts perpen- 
dicular to the longest common edge. Let us assume, as 
in Figure 2, that this edge lies in the y-direction. Then, 
we cut where the number of grid points covered by the 
intersections in each row changes. This list of rectangles 
is recursively tested for proper nesting. Obviously, this 
procedure is well suited to assign parents to each rect- 
angle at the same time. After having obtained a list of 
properly nested rectangles, they get information about 
their neighbors. 




FIG. 2. The result of the nesting procedure. 



FIG. 1. The effect of procedure saw up. 

To check, whether a rectangle is properly nested we 
calculate the sum of areas of intersections with all rectan- 
gles of the parent level. When this area equals the area of 
the actual rectangle it is guaranteed that this rectangle is 
properly nested. Otherwise, we proceed as follows. First, 
we determine the longest common edge of the just cal- 
culated intersections. Our strategy is to avoid coinciding 



2.5. Integral and local controls 

In order to extract physical properties of the simula- 
tion, it is necessary to calculate integral quantities like 
kinetic and magnetic energy as well as maxima of current 
density and vorticity. The latter are easily obtained by 
looping over all grids and all levels. Integral quantities 
are calculated in the following way. First, on the coarsest 
level the energy Ei eve i (swiss cheese energy) associated to 
the area not covered by grids of higher resolution is calcu- 
lated. This is repeated down to the lowest level. Finally, 
the energy is obtained as a sum over all energies Ei eve i . 
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2.6. Parallelization 

On shared memory machines our adaptive mesh refine- 
ment code can be parallelized in an effective and straight 
forward way. The main time of the program is spent in 
the procedure singlestep. Since the number of grids is 
much higher than the number of processors paralleliza- 
tion is done by distributing the grids to the processors. 
That means that as soon as a singlestep on a grid is fin- 
ished, the next grid is passed to the free processor. This 
results in a very effective utilization of all processors. All 
this can easily be done using standard Posix threads. The 
implementation on distributed memory machines using 
the shared memory access model is in work. 



3. NUMERICAL RESULTS 

In contrast to simulations of Frisch et al. |Q and 
Sulem et al. , we choose as initial condition a modified 
Orszag-Tang vortex, given by 

(p°(x, y) = cos{x + 1.4) + cos(y + 2.0) , 

iP a (x, y) = - [cos(2a; + 2.3) + cos(y + 6.2)] . 
o 

This initial condition, which was already used in turbu- 
lence simulations |l3|,|l(|, possesses less symmetry and is 
therefore more generic for the formation of small-scale 
structures. Computations are done with periodic bound- 
ary conditions on a square of length 2w. The initial spa- 
tial resolution was given by 256 2 grid points. 

The temporal evolution of the current density is shown 
in the contour plots of Figure |[ In addition to the con- 
tour levels, the rectangle hierarchy is plotted. The first 
plot shows the grid after the first refinement has taken 
place. The contour plot at time t = 2.2 contains already 
3 levels. At the final time t = 2.7 a total of 5 levels are 
present. Figure f| is a contour plot at the same time as 
the last one of Figure 0. To avoid hiding the sharpness of 
the current sheets no rectangles are included. In the ac- 
tual simulation, the refinement factor was equal to r = 2. 
On a workstation with 128 Mbyte of main memory, four 
refinements could be realized corresponding to a resolu- 
tion of 4096 2 grid points with a non-adaptive scheme. 
The limiting factor is the amount of main memory avail- 
able, whereas up to this resolution computational costs 
are very moderate. 

In the first picture the current sheets start to form, 
afterwards they evolve into thincr and thiner sheets and 
the maxima of current density and vorticity are increas- 
ing continually. The current density is growing expo- 
nentially in time. In Figure |B| a semilogarithmic plot of 
the maximum current density in the upper sheet is de- 
picted. Included is a fit to an exponential function given 
by jfit(i) = 0.5exp(2.115t). This functional behavior 
is in agreement with the results of Sulem et al. |ll[. A 
detailed analysis of the asymptotic scaling behavior and 





FIG. 3. 
and 2.7. 



Evolution of the current density at times 1.6, 2.2 
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FIG. 4. Current density at time 2.7. 



a comparison to the predictions in |T^] will be presented 
elsewhere. 
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FIG. 5. Amplification of current density. 



The second order upwind scheme produces substan- 
tial energy dissipation only if underresolved steep gradi- 



ents have formed. Therefore, the energy conservation is 
a measure whether the singular current sheets are suffi- 
ciently resolved. In Figure O we give a plot of energy as 
a function of time. To be more precise, total energy is 
conserved to within less than 1 %. 
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FIG. 6. Energy conservation. 
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In order to further illustrate that the current sheets are 
well resolved, in Figure |t] we show one dimensional cuts 
in cc-direction through the maximum of current density 
in the upper half of the integration range. In the upper 
plot the x-range equals the periodicity length. The lower 
one with a reduced plot range shows that the grid points 
of the finest levels very well resolve the current sheet. 
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FIG. 7. Cuts of current density in ^-direction through the 
maximum at time 2.72. 

In the previous section we mentioned that a refine- 
ment would take place when the discretization error for 
the nonlinearity exceeds a prescribed value e. The choice 
of the parameter e is crucial for the numerical accuracy. 
If e is taken too large, certain regions may be under- 
resolved which can lead to reconnection and violation of 
energy conservation. Decreasing systematically the value 
of e has the effect that reconnection phenomena are sup- 
pressed. Below a certain threshold the numerical results 
proved to be independent of e. The simulations shown in 
Figure |^ were performed with e = 0.025. 

Applying adaptive mesh refinement to the evolution 
of singular structures like current sheets in magnetohy- 
drodynamics is motivated by the expected reduction of 
memory needed to resolve them. This is well justified 
by the numerical results. To give the reader an impres- 
sion of how many grids are generated on the different 
levels and of the number of grid points contained in each 



level's grids, we display values for the level hierarchy at 
time t = 2.7 in the following tables. In Table I the results 
are shown for a simulation with refinement factor r = 2 
and in II for another one with r = 4. 



TABLE I 
Statistics for simulation with r = 2. 



level 


number of knots 


grid points 


1 T~l 1 PVpl 

111 1 v> V t/1 


n 


1 




70225 


l 


51 




168033 


2 


100 




341349 


3 


178 




734426 


4 


417 




1557221 




TABLE II 








Statistics for simulation with r = 4. 




level 


number of knots 


grid points 


in level 





1 




70225 


1 


49 




506073 


2 


195 


2331952 



From level to level the total number of grid points 
grows much less than by a factor of r 2 necessary for a 
non-adaptive treatment. For r chosen equal to 2, one 
can see that even for the very small value of e prescribed 
here it increases no more than by a factor of about 2. 
This promises that the compression rate will improve the 
more refinements are performed. 

In Table III simulations with different refinement fac- 
tors are compared with regard to the total number of grid 
points on all levels. The number of grid points on one 
data field with the same grid spacing as the finest level 
in the adaptive code is called the non-adaptive size. In 
the last row we give the ratio of the grid points, adap- 
tively and non-adaptively. For the investigated hierar- 
chy of 5 levels with refinement factor r = 2 this ratio is 
about 17 %. When the finest levels are equally resolved, 
the compression for both refinement factors is practically 
indistinguishable. For the comparison of adaptive ver- 
sus non- adaptive treatment, the compression rate based 
on counting grid points does not fully reflect the total 
improvement in main memory consumption. In upwind 
schemes several auxiliary fields have to be stored. In non- 
adaptive simulations these full sized fields are present all 
the time whereas here they are needed only temporarily 
during the execution of singlestep on a small grid. 



TABLE III 
Comparison of different refinement factors. 





r = 2 


r = 4 


total number of grid points 


2871254 


2908250 


non-adaptive size 


16851025 


16851025 


ratio 


0.170 


0.172 
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We want to finish this section with an impressive com- 
parison of the results for the amplification of the cur- 
rent density for several non-adaptive grid sizes and for 
the adaptive code. Figure || is a parametric plot of 
the maximum of current density as a function of the fit 
jtit(t) — 0.5exp(2.115i) already depicted in Figure [|. In 
addition to the results of the adaptive mesh refinement 
code we include data obtained with fixed grids of resolu- 
tions 128 2 , 256 2 and 512 2 . Until the simulations become 
underresolved, a linear behavior is also observed in the 
non-adaptive simulations. Then the upwind method in- 
troduces numerical viscosity leading to reconnection pro- 
cesses and substantial energy dissipation. 




20 40 60 80 100 120 



j flt (t) = 0.5exp(2.115t) 

FIG. 8. Comparison of adaptive and non-adaptive simula- 
tions. 

4. CONCLUSIONS 
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The complexity of adaptive mesh refinement compared 
to non-adaptive treatments should not be underesti- 
mated. On the other hand, the growing progress of ob- 
ject oriented programming languages helps enormously 
to reduce the difficulties in programming. To give some 
impression, the programs needed for regridding, nesting 
and the handling of data structures are only about 3000 
lines of C++ code. 

As we have demonstrated, adaptive mesh refinement is 
a powerful tool to study the evolution of singular struc- 
tures as the formation of current sheets in ideal MHD. 
Other problems of this type like in the axisymmetric |i~4| ] 
and the full three dimensional Eulcr equations are natu- 
ral candidates for this method. Work in this direction is 
in progress. 

Whether adaptive mesh refinement is also a useful con- 
cept for simulating turbulent hydro- and magnetohydro- 
dynamic flows will depend on how efficiently the small 
scale structures can be covered by hierarchically nested 
grids. 
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